home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / procssng / ccs / ccs-11tl.lha / lbl / hips / sources / getobj3d / connect.c next >
Encoding:
C/C++ Source or Header  |  1991-09-17  |  10.8 KB  |  487 lines
  1.  
  2. /* connect.c:   routines relating to the connectivity algorithm
  3.                 for getobj3d.c
  4.  
  5.     Brian Tierney,  LBL
  6. */
  7.  
  8. #include "getobj.h"
  9.  
  10. #define FASTER
  11. /* FASTER changes routines which only access the data
  12.    sequentially to use a faster method. The old method is retained
  13.    because it is a bit more readable. */
  14.  
  15. /***************************************************************/
  16. int
  17. identify_objects(thresh_value)    /* labels objects with a '1' */
  18.     int       thresh_value;
  19. {
  20.     int       num_on = 0, check_val;
  21.  
  22. #ifdef FASTER
  23.     register int i;
  24.     register u_char *bptr, *gptr;
  25.     register u_short *sptr;
  26.     register u_int *iptr;
  27.  
  28.     gptr = **grid;
  29.     if (pix_format == PFBYTE)
  30.     bptr = **image;
  31.     else if (pix_format == PFSHORT)
  32.     sptr = **sht_image;
  33.     else
  34.     iptr = **int_image;
  35.  
  36.     for (i = 0; i < nvoxels; i++) {
  37.     if (pix_format == PFBYTE)
  38.         check_val = (int) bptr[i];
  39.     else if (pix_format == PFSHORT)
  40.         check_val = (int) sptr[i];
  41.     else
  42.         check_val = iptr[i];
  43.  
  44.     if (check_val < thresh_value)
  45.         gptr[i] = 0;
  46.     else {
  47.         gptr[i] = 1;
  48.         num_on++;
  49.     }
  50.     }
  51.  
  52. #else                /* easier to read, about twice as slow  */
  53.     register int c, r, f;
  54.  
  55.     /* change all data to 0 or 1, depending on the thresh value */
  56.     for (f = 0; f < nframe; f++)
  57.     for (r = 0; r < nrow; r++)
  58.         for (c = 0; c < ncol; c++) {
  59.         if (pix_format == PFSHORT)
  60.             check_val = (int) sht_image[f][r][c];
  61.         else if (pix_format == PFINT)
  62.             check_val = int_image[f][r][c];
  63.         else
  64.             check_val = (int) image[f][r][c];
  65.  
  66.         if (check_val < thresh_value)
  67.             grid[f][r][c] = 0;
  68.         else {
  69.             grid[f][r][c] = 1;
  70.             num_on++;
  71.         }
  72.         }
  73.  
  74. #endif
  75.     return (num_on);
  76. }
  77.  
  78. /***************************************************************/
  79. void
  80. locate_surfaces()
  81. {                /* mark all surfaces as 2 */
  82.     register int c, r, f;
  83.  
  84.     /* label the surfacen points as 2 -- only checking 2-d slice */
  85.  
  86.     for (f = 0; f < nframe; f++)
  87.     for (r = 0; r < nrow; r++)
  88.         for (c = 0; c < ncol; c++)
  89.         if (grid[f][r][c] == 1) {
  90.             if ((r > 0 && grid[f][r - 1][c] == 0) ||
  91.             (r < nrow - 1 && grid[f][r + 1][c] == 0) ||
  92.             (c > 0 && grid[f][r][c - 1] == 0) ||
  93.             (c < ncol - 1 && grid[f][r][c + 1] == 0)) {
  94.             grid[f][r][c] = 2;
  95.             }
  96.         }
  97. }
  98.  
  99. /*********************************************************/
  100. void
  101. clear_background_image(bg_val)
  102.     int       bg_val;
  103.  /* sets the image pixels to zero based on the grid value */
  104.  
  105. {
  106.  
  107. #ifdef FASTER
  108.     register int i;
  109.     register u_char *bptr, *gptr;
  110.     register u_short *sptr;
  111.     register u_int *iptr;
  112.  
  113.     gptr = **grid;
  114.     if (pix_format == PFBYTE)
  115.     bptr = **image;
  116.     else if (pix_format == PFSHORT)
  117.     sptr = **sht_image;
  118.     else
  119.     iptr = **int_image;
  120.  
  121.     for (i = 0; i < nvoxels; i++) {
  122.     if (gptr[i] == 0) {
  123.         if (pix_format == PFBYTE)
  124.         bptr[i] = (u_char) bg_val;
  125.         else if (pix_format == PFSHORT)
  126.         sptr[i] = (u_short) bg_val;
  127.         else
  128.         iptr[i] = bg_val;
  129.     }
  130.     }
  131.  
  132. #else                /* easier to read, about twice as slow  */
  133.     register int c, r, f;
  134.  
  135.     for (f = 0; f < nframe; f++)
  136.     for (r = 0; r < nrow; r++)
  137.         for (c = 0; c < ncol; c++)
  138.         if (grid[f][r][c] == 0) {
  139.             if (pix_format == PFBYTE)
  140.             image[f][r][c] = (u_char) bg_val;
  141.             else if (pix_format == PFSHORT)
  142.             sht_image[f][r][c] = (u_short) bg_val;
  143.             else
  144.             int_image[f][r][c] = bg_val;
  145.         }
  146. #endif
  147.  
  148. }
  149.  
  150. /***************************************************************/
  151. void
  152. clear_background_grid()
  153.  /*
  154.   * at this point, all points in the desired object should be 3, so set any
  155.   * other points to 0
  156.   */
  157. {
  158.  
  159. #ifdef FASTER
  160.     register int i;
  161.     register u_char *gptr;
  162.  
  163.     gptr = **grid;
  164.  
  165.     for (i = 0; i < nvoxels; i++) {
  166.     if (gptr[i] == 3)
  167.         gptr[i] = 255;
  168.     else
  169.         gptr[i] = 0;
  170.     }
  171.  
  172. #else                /* easier to read, about twice as slow  */
  173.  
  174.     register int c, r, f;
  175.  
  176.     for (f = 0; f < nframe; f++)
  177.     for (r = 0; r < nrow; r++)
  178.         for (c = 0; c < ncol; c++) {
  179.         if (grid[f][r][c] == 3)
  180.             grid[f][r][c] = 255;
  181.         else
  182.             grid[f][r][c] = 0;
  183.         }
  184. #endif
  185. }
  186.  
  187. /***************************************************************/
  188. int
  189. mark_edge_neighbors(c, r, f, x_flag)    /* if neighbors are an edge, mark as
  190.                      * being part of object */
  191.     register int c, r, f, x_flag;
  192. {
  193.     int       cnt = 0;
  194.  
  195.     if ((r > 0) && (grid[f][r - 1][c] == 2)) {
  196.     grid[f][r - 1][c] = 3;
  197.     cnt++;
  198.     }
  199.     if ((r < nrow - 1) && (grid[f][r + 1][c] == 2)) {
  200.     grid[f][r + 1][c] = 3;
  201.     cnt++;
  202.     }
  203.     if ((c > 0) && (grid[f][r][c - 1] == 2)) {
  204.     grid[f][r][c - 1] = 3;
  205.     cnt++;
  206.     }
  207.     if ((c < ncol - 1) && (grid[f][r][c + 1] == 2)) {
  208.     grid[f][r][c + 1] = 3;
  209.     cnt++;
  210.     }
  211.     if (x_flag) {        /* also check x-slice */
  212.     if ((f > 0) && (grid[f - 1][r][c] == 2)) {
  213.         grid[f - 1][r][c] = 3;
  214.         cnt++;
  215.     }
  216.     if ((f < nframe - 1) && (grid[f + 1][r][c] == 2)) {
  217.         grid[f + 1][r][c] = 3;
  218.         cnt++;
  219.     }
  220.     }
  221.     return (cnt);
  222. }
  223.  
  224. /***************************************************************/
  225. int
  226. get_seed(sx, sy, sz)
  227.     int      *sx, *sy, *sz;
  228. /* looks for an object at the given seed point +-10 pixels:
  229.  * Starts at given seed, first look up and right 10 pixels, then
  230.  * looks down and left 10 pixels
  231. */
  232. {
  233.     register int c, r, f;
  234.  
  235.     int       bx, by, bz;
  236.  
  237.     fprintf(stderr, "\n Looking for object at (%d,%d,%d) \n", *sx, *sy, *sz);
  238.  
  239.     bx = *sx;
  240.     by = *sy;
  241.     bz = *sz;
  242.  
  243.     /* look at 4 edge neighbors to determine if good seed */
  244.     for (f = bz; f < bz + 10; f++)
  245.     for (r = by; r < by + 10; r++)
  246.         for (c = bx; c < bx + 10; c++)
  247.         if (grid[f][r][c] == 1) {
  248.             if (count_neighbors(c, r, f) > 3) {
  249.             *sx = c;
  250.             *sy = r;
  251.             *sz = f;
  252.             return (0);
  253.             }
  254.         }
  255.     /* try again */
  256.     for (f = bz - 10; f < bx; f++)
  257.     for (r = by - 10; r < by; r++)
  258.         for (c = bx - 10; c < bz; c++)
  259.         if (grid[f][r][c] == 1) {
  260.             if (count_neighbors(c, r, f) > 3) {
  261.             *sx = c;
  262.             *sy = r;
  263.             *sz = f;
  264.             return (0);
  265.             }
  266.         }
  267.     /* if get thru without returning */
  268.  
  269.     return (-1);
  270. }
  271.  
  272. /***************************************************************/
  273. int
  274. get_seed_guess(sx, sy, sz)
  275.     int      *sx, *sy, *sz;
  276.  /*
  277.   * if the user does not specify a seed value, then this routine located
  278.   * reasonable guesses for seeds. If the -a (find all objects) flag is
  279.   * specified, it begins the search from location (1,1,1). Otherwise it
  280.   * starts near the center of the data set, and if it doesn't find anything,
  281.   * then looks from location (1,1,1).
  282.   */
  283. {
  284. /* returns a seed for an object larger than 8 pixels */
  285.     register int c, r, f;
  286.  
  287.     static int bx = 0, by = 0, bz = 0;
  288.  
  289.     if (bx == 0 || by == 0 || bz == 0) {
  290.     if (find_all)
  291.         bx = by = bz = 1;    /* start at upper corner */
  292.     else {
  293.         bx = (ncol / 2) - 2;/* start near center of data set */
  294.         by = (nrow / 2) - 2;
  295.         bz = (nframe / 2);
  296.         if (bx < 1)
  297.         bx = 1;
  298.         if (by < 1)
  299.         by = 1;
  300.     }
  301.     }
  302.     if (!find_all) {
  303.     if (verbose) {
  304.         fprintf(stderr, "\n Looking for a seed value starting at: %d,%d,%d \n",
  305.             bx, by, bz);
  306.     }
  307.     for (f = bz; f < nframe - 1; f++)
  308.         for (r = by; r < nrow - 1; r++)
  309.         for (c = bx; c < ncol - 1; c++)
  310.             if (grid[f][r][c] == 1) {
  311.             if ((grid[f][r][c - 1] == 1) &&
  312.                 (grid[f][r][c + 1] == 1) &&
  313.                 (grid[f][r - 1][c] == 1) &&
  314.                 (grid[f][r + 1][c] == 1) &&
  315.                 (grid[f][r - 1][c - 1] == 1) &&
  316.                 (grid[f][r + 1][c + 1] == 1) &&
  317.                 (grid[f][r + 1][c - 1] == 1) &&
  318.                 (grid[f][r - 1][c + 1] == 1)) {
  319.                 *sx = c;
  320.                 *sy = r;
  321.                 *sz = f;
  322.                 return (0);
  323.             }
  324.             }
  325.     bx = by = bz = 1;    /* didn't find a seed, so try again from
  326.                  * upper corner */
  327.     }
  328.     if (verbose) {
  329.     fprintf(stderr, "\n Looking for a seed value starting at: %d,%d,%d \n",
  330.         bx, by, bz);
  331.     }
  332.     for (f = bz; f < nframe - 1; f++)
  333.     for (r = by; r < nrow - 1; r++)
  334.         for (c = bx; c < ncol - 1; c++)
  335.         if (grid[f][r][c] == 1) {
  336.             if ((grid[f][r][c - 1] == 1) &&
  337.             (grid[f][r - 1][c] == 1) &&
  338.             (grid[f][r][c + 1] == 1) &&
  339.             (grid[f][r + 1][c] == 1) &&
  340.             (grid[f][r - 1][c - 1] == 1) &&
  341.             (grid[f][r + 1][c + 1] == 1) &&
  342.             (grid[f][r + 1][c - 1] == 1) &&
  343.             (grid[f][r - 1][c + 1] == 1)) {
  344.             *sx = c;
  345.             *sy = r;
  346.             *sz = f;
  347.             bx = c + 1;
  348.             by = r;
  349.             bz = f;
  350.             return (0);
  351.             }
  352.         }
  353.     /* if get thru without returning */
  354.     if (verbose)
  355.     fprintf(stderr, " object not found \n");
  356.  
  357.     return (-1);
  358. }
  359.  
  360. /***************************************************************/
  361.  
  362. void
  363. init_grid(val)
  364.     int       val;
  365.  /* sets the entire grid to 'val' */
  366. {
  367.  
  368. #ifdef FASTER
  369.     register int i;
  370.     register u_char *gptr;
  371.  
  372.     gptr = **grid;
  373.  
  374.     for (i = 0; i < nvoxels; i++)
  375.     gptr[i] = val;
  376.  
  377. #else                /* easier to read, about twice as slow  */
  378.     register int c, r, f;
  379.  
  380.     for (f = 0; f < nframe; f++)
  381.     for (r = 0; r < nrow; r++)
  382.         for (c = 0; c < ncol; c++)
  383.         grid[f][r][c] = val;
  384. #endif
  385. }
  386.  
  387. /***************************************************************/
  388. int
  389. count_diag_neighbor_edges(c, r, f, x_flag)
  390.     register int c, r, f, x_flag;    /* max value returned is 12 */
  391. {
  392.     /*
  393.      * NOTE: this routine only counts neighbors that are marked '2'. There
  394.      * will also be some neighbors marked 0, but we don't want to count those
  395.      * because when labeling the edges we only look at the 2-D slice.
  396.      */
  397.  
  398.     int       neighbors = 0, tx, ty, tz;
  399.  
  400.     tx = ncol - 1;
  401.     ty = nrow - 1;
  402.     tz = nframe - 1;
  403.  
  404.     if ((r > 0) && (c > 0))
  405.     if (grid[f][r - 1][c - 1] == 2)
  406.         neighbors++;
  407.     if ((r > 0) && (c < tx))
  408.     if (grid[f][r - 1][c + 1] == 2)
  409.         neighbors++;
  410.     if ((r < ty) && (c > 0))
  411.     if (grid[f][r + 1][c - 1] == 2)
  412.         neighbors++;
  413.     if ((r < ty) && (c < tx))
  414.     if (grid[f][r + 1][c + 1] == 2)
  415.         neighbors++;
  416.  
  417.     if (x_flag) {        /* count x-slices too */
  418.     /* previous slice */
  419.     if ((f > 0) && (r > 0))
  420.         if (grid[f - 1][r - 1][c] == 2)
  421.         neighbors++;
  422.     if ((f > 0) && (c > 0))
  423.         if (grid[f - 1][r][c - 1] == 2)
  424.         neighbors++;
  425.     if ((f > 0) && (r < ty))
  426.         if (grid[f - 1][r + 1][c] == 2)
  427.         neighbors++;
  428.     if ((f > 0) && (c < tx))
  429.         if (grid[f - 1][r][c + 1] == 2)
  430.         neighbors++;
  431.  
  432.     /* next slice */
  433.     if ((f < tz) && (r > 0))
  434.         if ((grid[f + 1][r - 1][c] == 2))
  435.         neighbors++;
  436.     if ((f < tz) && (c > 0))
  437.         if ((grid[f + 1][r][c - 1] == 2))
  438.         neighbors++;
  439.     if ((f < tz) && (r < ty))
  440.         if ((grid[f + 1][r + 1][c] == 2))
  441.         neighbors++;
  442.     if ((f < tz) && (c < tx))
  443.         if ((grid[f + 1][r][c + 1] == 2))
  444.         neighbors++;
  445.     }
  446.     return (neighbors);
  447. }
  448.  
  449. /***************************************************************/
  450. int
  451. count_neighbors(c, r, f)
  452.     register int c, r, f;    /* max value returned is 6 */
  453. {
  454.     /*
  455.      * counts number of primary neighbors, used by 'get_seed'
  456.      */
  457.  
  458.     int       neighbors = 0, tx, ty, tz;
  459.  
  460.     tx = ncol - 1;
  461.     ty = nrow - 1;
  462.     tz = nframe - 1;
  463.  
  464.     if (f > 0)
  465.     if (grid[f - 1][r][c] == 1)
  466.         neighbors++;
  467.     if (f < tz)
  468.     if (grid[f + 1][r][c] == 1)
  469.         neighbors++;
  470.  
  471.     if (r > 0)
  472.     if (grid[f][r - 1][c] == 1)
  473.         neighbors++;
  474.     if (r < ty)
  475.     if (grid[f][r + 1][c] == 1)
  476.         neighbors++;
  477.  
  478.     if (c > 0)
  479.     if (grid[f][r][c - 1] == 1)
  480.         neighbors++;
  481.     if (c < tx)
  482.     if (grid[f][r][c + 1] == 1)
  483.         neighbors++;
  484.  
  485.     return (neighbors);
  486. }
  487.